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^ Abstract 

(N 

Approximate solutions of the Fisher equation obtained by different splitting meth- 
ods are investigated. The error of this nonlinear problem is analyzed. The order 

<^ of different splitting methods coupled with numerical methods of different order 

^ is calculated numerically and symbolically. 
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1 — 1 1. Introduction 

Splitting methods have been fruitfully used to solve large systems of partial 
differential equations. To find the exact solution of a given problem in practice 
is usually impossible. We can use numerical methods to obtain an approximate 
solution of the equations, although the discretized model can be still very diffi- 
i— h cult to solve. Reaction-diffusion models or models of transport processes have a 

structure that allows a natural decomposition of the equations, thus provide the 
opportunity to apply operator splitting schemes. Splitting methods help us reduce 
the complexity of the system and reduce computational time. With splitting it is 
possible to handle stiff terms separately and to solve each subproblem with a suit- 
able numerical method chosen to the corresponding operator. To solve a problem 
in practice we use operator splitting and numerical schemes which we will call the 
combined method. The use of operator splitting as well as the numerical methods 
result in some error in the solution. The error generated purely by splitting is 
called splitting error. This is the difference of the exact solution and the approx- 
imate solution obtained by splitting (assumed that we know the exact solutions 
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of the subproblems). Combined methods can generate both splitting error and nu- 
merical error. The study of this common effect on the solution is our main concern 
in this paper. Detailed study on the interaction of operator splitting and numerical 
schemes for linear problems can be found in [|Csom6s and Farago]. They clas- 



sify the errors that can occur using splitting methods and numerical schemes, give 
theoretical and numerical results on the order of the combined method for linear 
problems. Our aim is to characterize the error of this combined method therefore 
we calculate the order of the combined method for a nonlinear problem. We an- 
alyze the order of the error in the light of the characteristics of the splitting error 
and the numerical error. 

BSanz-Sernal and ULanser and VerwerTl discuss the splitting error in a general 
framework. The effect of operator splitting on the wave solutions of the Fisher 



equation is investigated by [Simpson and Landman |. Our aim here is to rigor 



ously analyze the interaction of splitting error and numerical error in the case of a 
nonlinear problem: the Fisher equation. The structure of our paper is as follows. 
In section 2 we introduce the basic idea of operator splitting in a general frame. 
In Section 3 we introduce the Fisher equation and recall some known results on 
it. We show how we apply splitting to solve the Fisher equation. In Section 4 
the splitting error is analyzed for reaction-diffusion problems in full generality. In 
Section 5 we calculate the order of the combined method for nonlinear problems 
in general. Section 6 contains the numerical results on the Fisher equation. 

2. Operator splitting 

Let us consider the following abstract Cauchy problem: 

U'{t)=A{U{t)) U(0)=U (1) 

with Uq £ X arbitrary. The set X is usually a space of functions with certain 
properties, U(t) £ X for every t ^ and A : X — >• X. Suppose that A can be written 
as the sum of two operators: A = A\ +A2. 

The most simple type of operator splitting is the sequential splitting. In this 
case the split problem is: 

U[(t)=A l (Ui(t)) £/i(0) = £/o (2) 
Ui(t)=A 2 {U 2 (t)) U 2 (0) = U l (t). (3) 
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The basic idea of splitting is to decompose the operator on the right hand side 
into the sum of simpler operators, and to solve the subproblems corresponding to 
the operators successively in each time step. More precisely, we solve the equa- 
tion only with operator A\ until time x (as if only the subprocess represented by A\ 
were present) and the solution in time x will be the initial condition of the equation 
with A2. It means that we return to the initial time and solve the equation with A2 
as well. The solution of the second equation in time x is called the approximate 
solution of the original problem in time x. This procedure is then repeated on the 
interval [t, 2t] etc. Thus, the simpler subproblems are connected to each other 
through the initial conditions. It is clear, that the numerical treatment of the sepa- 
rate subproblems is simpler. The most significant advantage of splitting is that we 
can exploit the special properties of the operators of the different subproblems and 
apply the most suitable numerical method for each of them. Thus we can obtain a 
more precise solution in a shorter time. 

We remark that the method can be used fruitfully in large models, for exam- 
ple global models of air pollution transport, or combustion or metabolic models, 
where the number of predicted variables is large and the number of the processes 
represented in the models is large. We refer to three works on air pollution models 
with application of operator splitting of [ |Lagzi et aL |. 

2.1. Splitting schemes 

We can define the different splitting methods by solving the subproblems suc- 
cessively in different orders and for different time lengths. The above described 
simplest scheme is called sequential splitting (SEQ). We solve the subproblems 
one after another using the same time length x, schematically S2(x)S{ (t), where St 
is the corresponding solution operator. In Marchuk-Strang splitting (MS) we usu- 
ally solve the subproblem with A 1 with time length x/2 then solve the other one 
with A2 for time length x and solve with A 1 again for time length x/2. Schemat- 
ically S{(x/2)S2(x)S\(x/2). We usually chose A2 to be the operator represent- 
ing chemical reactions. In general the operator that is stiff or nonlinear. In this 
given order we only need to solve the second subproblem once which can be 
of importance given the operator's properties. In weighted sequential splitting 
the solution in the next time step is a weighted average of the results of the two 
possible sequential splittings Si (t) £2(1) and S2(x)Si(x). In the special case of 
symmetrically weighted splitting (SW) we take the arithmetic mean of the results: 
(5i(t)5 2 (t) +S 2 (x)Si(x))/2. The extra work with MS and SW splittings bene- 
fits in second order accuracy compared to the first order of SEQ splitting. The 
nonsymmetric weighted splitting is of order one. In later sections we investigate 
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the SEQ, the MS and the SW splittings coupled with four different numerical 
methods, all of different orders. 



2.2. Splitting error, order of splitting 

We perform a semidiscretization on ([T]) in an equidistant manner with time step 
x. If we know the exact solutions of ([2]) and ([3]) we can generate an approximate 
solution to the original full problem ([j} in which error originated only from oper- 
ator splitting can arise. If we denote the exact solution by U and the approximate 
solution by U then the local error of operator splitting is 

E(x) :=U(x)-U(x) 

Both solutions start from the common initial value and after time x the difference 
E(x) is called splitting error. Naturally the splitting error can be defined at any 
point of time during integration, if U(t) = U (t) then U (t + x) — U (t + x) is the 
local error at t. 

For linear operators it is easy to show by Taylor expansion that the local order 
of SEQ equals 1 since the error becomes E(x) = Kx 2 + 0(x 3 ). For nonlinear op- 
erators we need the definition of the Lie-operator and we can perform the analysis 
with Taylor expansion using the Lie-operators. We refer to ULanser and Verwerll 
for detailed derivation of the nonlinear case. From the literature on operator split- 
ting it is well known that the MS provides second order accuracy, so does the SW 
splitting, see [Farago and HavasiJ . 



2.3. Splitting of reaction-diffusion equations 

In the case of reaction-diffusion equations there is a natural decomposition 
of ([T]). The operator A\ represents the process of diffusion and A2 the chemical 
reactions. In this case A\ is a linear and unbounded operator, A2 is usually a 
nonlinear operator. If we have M species and N denotes the spatial dimension of 
the problem then U £ C l (R + ,C 2 (R N ,R M )) that is for a given time point t £ K+ 
the function U maps the concentration of all species for every given point in space 
M. N , so it is a function of x £ M. N . In other words U : t (ui(t,x),...,UM(t,x)), 
where w;(V,x) is the concentration of the zth species that is the spatial distribution 
of the zth species. Then the function spaceX := C 1 (M + ,C 2 (M W ,M M )). 

3. The Fisher equation 

The Fisher equation is: 

J d t u(t,x) = d 2 u(t,x) + u(t,x)(l — u(t,x)) x£lR,?>0 



u(0,x) = rj(x). 



(4) 
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There is only one chemical species present and one spatial variable here. This 
equation was originally derived to describe the propagation of a gene in a popu- 
lation [Fisher]. It is one of the simplest nonlinear models for reaction-diffusion 
equations. Such equations occur, e.g., in combustion, mass transfer, crystalliza- 
tion, plasma physics, and in gen eral phase tran sition problems. See a discussion 



on reaction-diffusion models in IIErdi and Tothft and [ Murray) . For the initial con- 
dition: 

u(0,x) 



(l+kexp(x/V6)) 2 
wave form solution of the equation is known: 

1 

u(t,x) 



(l+kexp(-^t+^y/6. 



x, 



and for: 



u(0,x) 



1 



u(t.x) 



(1 +kexp(-x/V6)) : 
1 



(l + ^exp(-|t-i v / 6x))2 

We investigate three different splitting methods applied in the solution of this 
equation. A natural way to split the Fisher equation is to decompose it into two 
subproblems: one for the diffusion and one that corresponds to the reaction part 
of the right hand side. Thus the definitions of the subproblems are: 

d t m(t,x) = d*ui(t,x) 

Ul(0,*) =T?lW 
d t U2(t,x) = U2(t,x)(l—U2(t,x)) , 

u 2 (0,x) = tj 2 (x), 

where the initial condition 7]2( x ) — u\(l,x) connects the equations. The lower 
indexes help distinguish the solutions of the different problems. We define the 
operators A\ and A 2 as follows: A\u := d^u, A2(u) := w(l — u). These operators 
are independent of time. Following the convention for linear operators we neglect 
the parenthesis in Aim, at the same time we would like to emphasize that A2 is 
nonlinear: A 2 (u). 
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Figure 1: The exact solution of (6l, rj2{x) — — sin(x) + 1, t € [0,1] and x € [0,4tt]. The same 
initial condition is used as in section 6. 



The exact solution of problem (|6]) is known, it is: 



Since 

1 when 772 (jc) t^O 



lim »2(/..v) — , „ . , „ 

/v ' y [0 when r7 2 (x) = 

the solution has two stationary states, namely: U2(t,x) = 0, U2(t,x) = 1. The 
U2(t,x) = 1 solution is asymptotically stable, whereas zero is an unstable equilib- 
rium. 

Knowing the exact solution of this subproblem as a function of the initial 
condition means that we can symbolically solve this subproblem in each time 
step during the splitting procedure. It might be worth using the exact solution for 
comparisons in the study of the effect of splitting. The exact solution of ([5]) is of 
no real use. 



4. Commutation of diffusion and reaction 

In this section we investigate the conditions under which the SEQ has zero 
splitting error for reaction-diffusion systems in general. In reaction-diffusion 
equations there are two operators present on the right hand side, a linear and a 
nonlinear operator. Let us consider 
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U'(t)=A{U(t))+R(U(t)) U(0) = U , (8) 

where U : Rq -> X with X = {<p = (<pi,..., 9m); <Pm : R 3 -> Rq > m = 1 5 — >^}- 
We have M different species undergoing the processes of diffusion and reactions. 
The operator A represents the diffusion assuming that there is no cross diffusion 
present. In the above vectorial form A is an operator matrix essentially with the 
spatial Laplacian in the diagonal and zeros everywhere else. A is a linear operator. 
The operator R usually acts as compositions with the multivariable polynomials 
Ri,R2--.,Rm, it is a nonlinear operator. That is 



R(<p) = (R l ((p h ...,(p M ),...,R M {(p h ...,(p M )) 



and 



A((p) =A((p h ...,(p M ) = (DiA(p u ...,D M A(p M ). 
The global equation ([8]) in the local form is: 

U'(t)(x) =A(U(t))(x) +R(U(t))(x) U(0)(x) = U (x) 

Therefore, another form of ([8]) is 

d t m(t,x) = DiA«i(f,x)+i?i(u(f,x)) 

d t u M (t,x) = D M Au M (t,x)+R M (n(t,x)) 



(9) 



with D m > 0, m = 1, . . . ,M, u(-) = (wi(-)j ... : um(-)) and x G M 3 . Since A is linear 
its derivative is A itself. The derivative of R := (R\,...,Rm) which is the Jacobi 
matrix: 

/ diRi ... d M Ri \ 
R'= : : 

y d\R M ... 3 m Rm j 

In the following we derive the condition of zero splitting error. The sufficient 
condition of zero splitting error is [Hundsdo rfer and VerwerB that for every <p e X 



(A'oR-R'oA) (<p) =0 

Upon applying this we get 

M 

= D m AR m ((p(t,x)) - £ d k R m ((p(t,x))D m A(p k (t,x) 



(10) 
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For the sake of simplicity we will not carry the argument (f,x), it will no lead to 
misconceptions. 

M M M 

= d x £ d k R m ((p)d x (p k + d y £ d k R m ((p)d y (p k + d z £ d k R m ((p)d z (p k - 

k=\ k=\ k=\ 

M 

-Y, d k R m(<p)&<Pk m=l,...,M. 

Each of the first three terms is a sum of products thus, according to rule of differ- 
entiation of a product one gets: 

M M M 

£ d x (d k R m ((p))d x (p k + £ dy(d k R m ((p))d y (p k + £ d z (d k R m (<p))d z (p k + 
k=\ k=\ k=i 

M M M M 

+ £ d k R m ((p)d^(p k + £ d k R m ((p)d^(p k + £ d k R m ((p)d?(p k - £ d k R m ((p)A(p k 
k=i k=i k=i k=i 

Since the 4th, 5th, 6th and 7th terms eliminate each other, only the first three terms 
remain. Upon performing the differentiation the above expression becomes: 

MM MM 

£ £ djd k R m ((p)d x (pjd x (p k +Y, E djd k R m ((p)d y (pjd y (p k + 

k=\ j=l 7=1 
M M 

+ £ £ djd k R m ((p)d z (pjd z (p k 
k=lj=l 

Using the notation: d x (p := {d x <p\, , <9 x <Pm) we can reformulate this as: 



= < ^((P) • d X (p, d x (p> + < <(</)) • ^(P, dy? > + < ^((P) • ^ Z (P, drf > . 

The splitting error is zero if the above expression equals zero. Using the notation: 

/ d x (pi d y (pi <9 z <pi 



d(p:=(d x (p,d y (p,d z (p) = 



\ d x (p M d y <pM d z (p M 



we can write the formula above in the short form: 



(R'M-dcp) T .dcp = 0. 
Now we can formulate a statement. 

Theorem 1. With the notations above the error ofSEQ is zero if for every function 
<p and for all m= 1,2, ....M (/?^(<p) • d(p) T -d(p = holds. 

Remark 1. We get the necessary condition of zero splitting error if we require 
condition ( [70] ) to hold only for the solution function of ([§]). Naturally without 
knowing the exact solution we can not check whether this holds or not. But the 
above formula provides us a sufficient condition for zero splitting error. If the 
equation holds for every possible function (p then it will hold for the solution as 
well. This holds for every function iff all the entries ofB!' m is zero, which means that 
R m is a polynomial of at most first degree, m= 1, ...,M : we only have first order 
reactions. Most of the practical problems have reaction terms of higher order, 
therefore there is almost always a splitting error. Our aim here is to examine the 
effect of splitting error in combined methods. 



Remark 2. Condition (10) is sufficient in the case of MS and SW splittings as 



well. Without going into the details condition (flO]) ensures that e TA e TR = e T ( A + R 



where A and R are the Lie-operators of A and R. In the case of MS we need 
e2 A e rR e2 A = g T ( A+R ) which obviously holds if e tA e rR = <? T ( A+R ). Theorem 1 
remains valid in the case of MS and SW splittings. 



5. Order of combined methods 

When we solve partial differential equations we can use some kind of split- 
ting but we can not avoid applying some numerical method as well. So in prac- 
tice we use a combined method, a mixture of operator splitting and a numerical 
scheme and generate a solution for a nonlinear partial differential equation like 
Q. We use the Taylor-formula to determine the order of the local error of this 
combined method. The Taylor-formula in normed vector spaces can be found in 
e. g. [KomornikJ. Here we recall the Taylor-formula for normed vector spaces: 

Theorem 2. Iff : X — > Y is n times differentiable in a EX and h — >■ 0, then 
f(a + h) = t^P-h k + e(h)\\h\\ n 

where h k := (h,...,h) G X k and lim e(h) = 0. 
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Suppose that U is the solution of the equation: 



U'(t)=AU(t)+R(U(t)) (11) 

From now on in this section we restrict our investigations to problems where 
A : X X is a bounded linear operator defined on the whole set X, therefore 
the following derivations are not directly applicable to the case of diffusion. A is 
differentiable and its derivative is A itself for every x E X. The operator R.X^X 
will act as a composition with a differentiable nonlinear function thus R is a differ- 
entiable mapping as well. Based on ( fTTj ) and the chain-rule U' is a differentiable 
function and U'(t) E X. 

The Taylor-expansion of U in time to is 

£/ T := t/(* + t) = U(t ) + U'(t )x+ l -U"(to)% 2 + e(r) || t|| 2 . 

The norm we will neglect from now on since U : R — > X, T denotes a positive real 
number. U'(to) is given by ( fTT| ), we get U"(to) by differentiation of ( fTT| ): 

t/"(*o) =A , (t/(?o))ot/ , (?o)+ J R , (t/(fo))t/ , (?o) =A(t/ , (?o))+i? , (t/(?o))t/ , ao) = 

= A(A([/(* )) +i?(f/(?o))) +i? , (t/(?o))A(f/(? )) +i? / (C/(f ))^(f/(fo))- 

U (to) is denoted by Uq. Using this notation if the value of U (to) = Uq is known 
then we can approximate U t : 

U t = Uo+(A(Uo)+R(U ))t+ (12) 

+ l - (a(A(U )) +A(R(U )) +R'(U )A(Uo) +R , (Uo)R(Uo)^j x 2 + c(t)t 2 . 

Beyond the conditions we already mentioned this form of U z exists if Uq E D(A 2 ) 
and R(Uq) E D(A) which naturally hold if A : X -»■ X and 7? : X ->■ X. For the the- 
orems in this section we will need 7? to be three times continuously differentiable. 

5.1. Methods of first order 

Theorem 3. The sequential splitting combined with the first order Euler forward 
scheme provides a first order method. 

The proof will be given in two steps. 
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5.1.1. Linear-Nonlinear 

If we use SEQ starting with the nonlinear problem corresponding to R com- 
bined with Euler forward method for both subproblems we get: 

f V(x) = U + xR(U ) 
\ U t = V(x) + xA(V(x)) 

Proof. The approximation of the solution in time x is: 

Z = V(x) + xA(V(x)) = U + xR(U ) + xA(U + xR(U )). 
Since A is linear we have: 

U T = U + xR(U ) + xA(U ) + x 2 A(R(U )). (13) 
The local error generated in this step of length x based on { fT2] ) and ( fT3| ) is: 

2 

U Z -0 Z = (a(A(U )) -A(R(U )) +R'(U )A(Uq) +R'(Uo)R(U )^ % - + e(x)x 2 . 

□ 

5.1.2. Nonlinear-Linear 

If we use SEQ starting with the linear problem corresponding to A combined 
with Euler forward method for both subproblems we get: 

f V(x)=U + xA(Uq) 

1 U z = V(x) + xR(V(x)). 

Proof. The approximation of the solution in time x is: 

U z = V(x) + xR(V(x)) = U + xA{Uq) + xR(U + xA(U )). 

Let us define the function F : R — > X in the following way: F(8) := R(Uq + 
8xA(U )). Then 

F(0) =R(U ) andF(l) = j R([/ + tA([/ )), 

and F is differentiable since it equals to R o f with :=Uq + 8xA(Uq) (where 
is differentiable), and according to the chain-rule 

F'{8) = R'(U + 8xA(Uo))xA(U ) F"(8) = R"(U + 8xA(U ))(xA(U )) 2 

and 

F (")(8)=R^(Uo + 8xA(Uo))(xA(Uo)) n 
For the Taylor-expansion of F we need a similar but more specific relation. 
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Proposition 1. IfF : R — >■ X is n times differentiable in every point of [a,b] then 
there exists c € (a,b) such that 

In other words 

with ||£i (c) || < 1. The expansion of F around gives: 

F(l)=F(0)+F'(0) + E l (c) l -\\F"(c)\\ ce [0,1] 

implying 

J?(tfo+TA(tf ))=*(tfo)+#(tfo)TA(^ 

= tf(C/ ) +^(C/ )tA(C/ ) + 6(t)t, 

where 

e(x) = e l (c)^\\R"(U + cxA(U ))A(U ) 2 x\\. 

This tends to zero if x tends to zero. This is the only relevant property of £ here. 
Although c can change as x changes but since ||£i(c)|| ^ 1 we can ignore e's 
dependence on c through £i . 



U T = U + xA(U ) + x (R(U ) +R'(U )xA(U ) + £(x)\\x\ 

= U + xA(U ) + xR(U ) +R\Uo)x 2 A(U )+£(x)x 2 , 
U T = U + x(A(U )+R(U ))+ 

2 

+ (a 2 (U ) +A(R(U )) +R'(Uo)A(U ) +R'(Uo)R(U )) % -+e(x)x 2 . 
Comparing this with the approximation we get: 

2 

U T -U T = (a 2 (U ) +A(R(U )) -R'(U )A(Uq) + R'(U )R(U )) % - +e(x)x 2 

□ 
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5.1.3. Weighted Splitting 

Theorem 4. SW combined with the first order Euler forward method provides a 
method of first order. 

Proof. Here we simply use the above results with some CO G [0, 1] parameter: 
U t = co(U + t(R(U )+A(U )) + t 2 A(R(Uo))) + 

HI-oj)(Uo + t(A(Uo)+R(Uo))+R'(Uo)t 2 A(Uo)+R''(Uo)t 3 (A(U )) 2 ^ + £(t)t 2 , 

U T = Uo + t(A(U ) +R(U )) + (coA(R(U Q )) + (1 - to) (R f (U )A(U ))) t 2 + c(t)t 2 . 
For the local error we have: 

U T -U Z = 

= (a 2 (U ) + (1 - 2<o)A(U )R(U ) + (2w - \)R'(U )A(U ) + R' (U )R(U )^ y + 

+£(T)T 2 , 

for CO = - we have 
2 

U r - U T = (a 2 (U ) +R'(Uo)R(U )) y + £(T) t 2 . 

□ 

Conclusion: Although the SW is of second order its combination with the first 
order Euler method provides only first order accuracy. 

5.2. Methods of higher order 

The above derivation can be used to determine the order of combined methods 
with higher order numerical schemes. The method can be extended for schemes 
of arbitrary order although the calculations become very complicated as the order 
increases. As an example let us consider the improved Euler scheme which is of 
second order and combine it with SEQ: 

Theorem 5. The second order improved Euler scheme combined with SEQ pro- 
vides a first order method. 

Again, the proof will be given in two steps. 
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5.2.1. Nonlinear-Linear 

V(x)=U + xA(U + ^A(U )) 
U, = V(x) + xR(V(x) + ^R(V(x))) 
Proof. The approximation of the solution is in time x is: 

U x = V(x) + xr(v(x) + !*(V(t))) = U + xA(U + X -A(U )) + 

+xR (U + xA (U + X -A(U Q )) + X -R (u + xA (U + X -A(U ))) ) = 

x 2 

= U + xA(U ) + -A 2 (U )+ 

+xr(u + x(a(U ) + X -A 2 {U Q ) + l -R{U Q + xA(U ) + yA 2 (t/ )))) = 

The underlined part is the coefficient of 8 in the argument of F. The first order 
Taylor-expansion gives: 

x 2 

= U + xA(U ) + -A 2 (U ) + xR(U Q ) + 

+xR'(U q )x(a(Uq) + X -A 2 (U Q ) + ^R(U + xA(U ) + yA 2 (t/ ))) +e(x)x 2 = 

= U + x(A(U )+R(U )) + 
+ y (A 2 (Uo) + 2R'(U )A(U ) +R f (U )R(U + xA(U ) + yA 2 (t/ ))) + 

T 3 

+-R\U )A 2 (Uq) + £(x)x 2 = 
Taking the Taylor-expansion again, the coefficient of 8 is the underlined part: 

= U + x(A(U )+R(U )) + 

2 2 

+ X - (a 2 (Uq) +2R'(Uo)A(U ) +R'(U ) (R(Uo) + xR'(Uo)A(U ) + X ^R'(Uo)A 2 (U )))) + 

+e(x)x 2 = 

14 



= Uo + t(A(Uo)+R(Uo)) + ^-{a 2 (Uq)+2R\Uo)A(Uo)+R'(U )R(U ))+£(t)t 2 
U T -Ur=(A(R(Uo))-R'(Uo)A(Uo)^+£(t)T 2 . 

□ 

5.2.2. Linear-Nonlinear 

The proof of the linear-nonlinear case is more straightforward: 
Proof. 

V(x) = U + xR(U +^R(U )) 
U z = V(x) + xA(V(x) + ^A(V(x))) 

U t = V(x) + xa(v(x) + ^A(V(x))) = 

= U + xR(U +^R(U )) + 

+xa(u + xR(U + \r(U q )) + % -a(u q + xR(U + ^o)))) = 

= U + x(R(U ) + ^R(U )R'(U Q )) + xA{U )+ 

+T 2 (A(7?(t/ ) + |^(W(t/o)) + \a 2 (U )) +£(t)t 2 = 

2 

= U + x(R(U ) +A(U )) + ^(r(U )R'(Uo) + 2A(R(U )) +A 2 (U Q )^+e(x)x 2 

2 

U x -U x = (r'(U )A(U q ) -A(R(U ))) y + £(t)t 2 . 

□ 

As we can see this combined method is of first order. Although the applied 
numerical scheme ensures second order accuracy the use of sequential splitting 
results in order reduction. We followed the same ideas and calculated the orders 
for combinations of the introduced splitting methods and four different numerical 
schemes. The table below contains our results on orders of different splittings 
coupled with different numerical methods. Symbolic calculations on for example 
MS splitting coupled with 4th order Runge-Kutta method becomes complicated. 
An algorithm was written in Mathematica for these symbolic calculations. 
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exp. Euler (1) 


impr. Euler (2) 


Heun (3) 


Runge-Kutta (4) 


SEQ (1) 


1 


1 


1 


1 


SW (2) 


1 


2 


2 


2 


MS (2) 


1 


2 


2 


2 



Table 1 : Local orders of combined methods for ( fTT| 



The order of the methods are in the parenthesis. A study of the order of com- 
bined methods for bounded linear problems can be found in [ |Csom6s and Fara gd | . 
Their results say that the order of the combined method is s := min{p,r}, if p is 
the order of splitting and r denotes the order of the numerical method. The num- 
bers of the above table are in accordance with their results. 



6. Numerical experiments 

Here we introduce our numerical results on the Fisher equation. We solved 
both subproblems ([5]), ([6]) using numerical methods of four different orders; the 
explicit Euler, the improved Euler method which is of second order, the third 
order Heun and the fourth order Runge-Kutta method. We investigate the SEQ, 
SW and MS splitting methods which are of first and second orders. We calculated 
the errors and orders of these combined methods numerically. Our test problem is 
the following initial-boundary value problem: 

d t u(t,x) - 
u(0,x) - 
u(t,0) - 
u(t, An) - 

where x e [0,4n] and t e [0, 1] 



d^u(t,x) + u(t,x)(l 

l+0.9sin(» 

1 

1 



u(t,x)) 



(14) 



We performed a spatial semidiscretization with length parameter Ax = |^ that 
is we divided [0, An] into N = 30 parts of equal length. Our tests showed that 
finer divisions provides no significantly more accurate solutions that is the ob- 
tained error is of the same magnitude as with N = 30. We approximated the spatial 
derivative with the well known second order scheme: 

2 , ^ u(t,X i+ i)-2u(t,Xi) + u(t,Xi-i)) 
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Figure 2: Reference solution generated by fourth order Runge-Kutta scheme, T = 0.01. 



After temporal discretization with x = 0.01 we solved the full problem ( |T4| ) with 
the fourth order Runge-Kutta method. Taking a smaller time step resulted in 
solution that differs only in a magnitude of 10 6 . This provides the reference 
solution for our study. In the experiments we used spatial division N = 30 in 
every case, in fact we investigated the convergence of the semidiscrete submodels 
to the semidiscrete model: the reference solution. On the connection between the 
convergence to a semidiscrete model and convergence to a continuous model see 
ULarsson and Thomeell . 

6.1. Determination of the local order 

The local error E{x) of a method is of order s if 

E(t) = 0(t s+1 ) 

or by using a different formulation the order is s if 

E(T) 



lim 



V 

is finite where s G N is the smallest number with this property. Thus in practice 
we can estimate the order of the combined method by calculating this limit with 
different fixed values of s until we find the appropriate one. For a fixed numerical 
step size h the limit means naturally that T — > /z, therefore we choose h to be as 
small as possible close to the smallest number that our computer can represent. 
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Another way to calculate the order of the method is the following. Considering 
the two formula above we can conclude that 

E{x) 

T ,+ l ~ L 

for small T-s, where c is a constant which does not depend on t. We can take the 
logarithm of both sides: 

log£(t) « (s + l)logr + logc. 

This defines a straight line whose steepness gives the order of the method. The 
table below contains results for the local order s. The order of splitting is p and 
r is the order of the numerical scheme. The time step of the numerical method h 
was the tenth of splitting time step T. 



s {h = 0.H) 


exp. Euler (r = 1) 


Heun (r = 3) 


Runge-Kutta (r = 4) 


SEQ (p=l) 


0.98 


0.98 


0.98 


SW (p = 2) 


0.83 


1.99 


1.99 


MS (p = 2) 


0.93 


1.96 


1.96 



Table 2: Orders of combined methods for ( fl4] i 



It is interesting that with the explicit Euler method increasing the order of 
splitting does not improve the results. 

6.2. Estimation of the global order 

We used the time steps T ; in such a way that the evaluation of the error is 
simple, since the corresponding division of the time interval is a subset of the one 
of T = 0.01. 





T 2 


^3 


T 4 


T 5 


T 6 




0.2 


0.1 


0.0625 


0.05 


0.04 


0.025 


0.02 



As for the error we know that: 

£(Ti)«C-Tf 
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for small T-s, where c is a constant which does not depend on x. So we can write: 



E{% 2 ) \x 2 
We can take the logarithm of both sides: 

For each i: 



log— r ^plog 



E(% i+ \) 

So 

log 

log^ ~ P ' 

Evaluation of the left side shall give us the same value for every i= 1 , 2, 3, 4, 5, 6. 
The following table contains the results of this calculation for different splittings 
and numerical methods. 



p(h = T) 


r=l 


r = 2 


r = 3 


r = 4 


SEQ (p=l) 


1.04 


0.99 


1.08 


1.08 


SW (p = 2) 


1.02 


2.07 


2.01 


1.98 


MS (p = 2) 


1.02 


2.07 


1.95 


1.998 



Table 3: Orders of combined methods for ([14) 

Here we had the same experience as with the calculations shown in Table 2. 
Since the solution of ([6]) is given in |7]) as a function of the initial condition we can 
use it in calculations instead of the numerical solution. The table below contains 
results generated by using ([7]) in each time step. 

The first column fits into the general scheme. 
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p(h = x) 


r=l 


r = 2 


r = 3 


r = 4 


SEQ (p=l) 


1.03 


1.02 


1.01 


1.01 


SW (p = 2) 


1.01 


2.06 


1.95 


1.98 


MS (p = 2) 


1.03. 


2.00 


1.99 


1.99 



Table 4: Orders of combined methods for ( 14 1, using dTJi 



6.3. Splitting into three operators 

Above we saw that the splitting error is zero if the reaction term is of first 
order. Another, natural decomposition of the right hand side of the Fisher equation 
is when we separate the reaction part into two terms. The split problem is the 
following: 

J d t u\{t,x) = d^ui(t,x) J d t U2(t,x) = U2(t,x) J d t u 3 {t,x) = —u^(t,x) 
\ki(0,x) =rj 1 (x) \u 2 (0,x) =7] 2 (x) \u 3 (0,x) = tj 3 (jc) 

(15) 

Then A\u = d^u, A 2 u = u, A 3 (u) = —u 2 . Since A 2 is a first order polynomial the 
splitting error of the SEQ S1S2 is zero according to theorem 1. It is easy to prove, 
that it is not zero for 53 coupled with any of the other two operators. 



s (h = t) 


r = 1 


r = 2 


r = 3 


r = 4 


SEQ (p=l) 


0.96 


1.96 


2.97 


3.96 


SW (p = 2) 


0.96 


1.96 


2.97 


3.96 


MS (p = 2) 


0.96 


1.96 


2.96 


3.96 



Table 5: Orders for splitting of d,u — d^u + u 



In the solution of subproblems associated with operator A 2 and A3 we used the 
exact solution. The first subproblem was solved numerically. Considering the MS 
type splittings S3S 2 SiS 2 S3 and S 2 S3S\St,S 2 it is reasonable to expect more accu- 
rate solutions with in the case St,S 2 S{S 2 S 3 since the neighbors Si and S2 generate 
no splitting error. Whereas S 2 Si,SiS?,S 2 generate splitting errors between every 
neighboring operators. Figure 5 and figure 6 proves that although both MS pro- 
vides the expected second order accuracy (combined with a third order or fourth 
order numerical schemes) the hypothetic relation in accuracy turns out to be the 
opposite. 
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(T,+r, +l )/2 



Figure 3: Order approximation of the S1S2S3 (squares) versus S1S3S2 (discs) SEQ type 
splittings. Generated by the fourth order Runge-Kutta scheme with time steps T = 
0.1, 0.0625, 0.05, 0.04, 0.025, 0.02. 



1.2 1.3 1.4 1.5 1.6 1.7 



Figure 4: Order approximation of SEQ type splittings as the function of t,/t, 



i+l- 
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order 
2.00 r 



(T,+r, +l )/2 



Figure 5: Order approximation of the SW splitting which is the arithmetic mean of the six possible 
sequential splittings. Generated by the fourth order Runge-Kutta scheme. 



order 
2.00 r „ 



(Ti+r, + i)/2 



Figure 6: Order approximation of the S3S2S1S2S3 (discs) versus S2S3S1S3S2 (squares) MS type 
splittings. Generated by the third order Heun scheme. 



7. Discussion and perspectives 

We presented symbolic calculations for orders of PDE solving methods. Our 
motivation is to predict the order in the case when beside numerical procedures of 
certain order operator splitting is also used. We calculated the order of combined 



methods applied for nonlinear PDE-s like ( |TT| ), where a bounded linear operator 
and a nonlinear operator is present. We presented numerical calculations on a 
test problem with the diffusion operator which is an unbounded linear operator. 
Although the results are in accordance with our theoretical results the methods 
used in section 5 are strictly correct for the case of bounded linear operator. The 
results of section 6 indicates that the combined method inherits the smaller one of 
the order of the splitting and the numerical method, the extension of the methods 
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order 
2.00 - 



1.2 1.3 1.4 1.5 1.6 



Figure 7: Order approximation of MS type splittings as the function of T//t,-+i. 



(t,+t m )/2 



Figure 8: Order approximation of the 5352 Si 5253 (discs) versus 52535i5352 (squares) MS type 
splittings. Generated by the fourth order Runge-Kutta scheme. 



1.! 



1.1 1.2 1.3 1.4 1.5 1.6 1.7 



Figure 9: Order approximation of MS type splittings as the function of t,/t, 



+ i ■ 
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used in section 5 to unbounded operators is not obvious. Our main focus is on 
reaction-diffusion equations so we plan to find a method which allows us to repeat 
the results of this paper for unbounded linear operators. 

We also intend to extend the results of | Csomos and Farago | to nonlinear prob- 
lems. Besides we are working on reaction-diffusion simulations on the sphere. In 
the future we are going to apply the presented methods to practical problems e. g. 
in the simulation of combustion. 
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